Multiple extinction routes in stochastic population models 
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Isolated populations ultimately go extinct because of the intrinsic noise of elementary processes. 
In multi-population systems extinction of a population may occur via more than one route. We 
investigate this generic situation in a simple predator-prey (or infected-susceptible) model. The 
predator and prey populations may coexist for a long time but ultimately both go extinct. In the first 
extinction route the predators go extinct first, whereas the prey thrive for a long time and then also 
go extinct. In the second route the prey go extinct first causing a rapid extinction of the predators. 
Assuming large sub-population sizes in the coexistence state, we compare the probabilities of each of 
the two extinction routes and predict the most likely path of the sub-populations to extinction. We 
also suggest an effective three-state master equation for the probabilities to observe the coexistence 
state, the predator-free state and the empty state. 
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I. INTRODUCTION 



Isolated populations ultimately go extinct with prob- 
ability one, even in the absence of adverse environmen- 
tal variations, because of the presence of an absorbing 
state at zero population size and the intrinsic noise of el- 
ementary processes. Population extinction risk is an im- 
portant negative factor in viability of small populations 
U [1[ , whereas extinction of an endemic disease from a 
community [l| is of course a good thing. For decades, 
quantitative analysis of population extinction, caused by 
intrinsic and extrinsic noises, has been in the focus of at- 
tention from bio-mathematicians, see e.g. Ref. p] for a 
review. More recently, it has also attracted much atten- 
tion from physicists |3l-|24|. as a highly relevant example 
of a rare large fluctuation far from thermal equilibrium. 

In multi-population systems extinction of a population 
may occur via more than one route. In this paper we 
analyze this generic situation in a simple predator-prey 
model where the predators need the prey for survival. 
The same model can describe spread of an infectious dis- 
ease in an isolated community. Although the predator 
and prey populations (or infected and susceptible pop- 
ulations) may coexist for a long time, they ultimately 
both undergo extinction, and this happens via one of two 
routes. In the first route, that can be called sequential, 
the predators go extinct first, whereas the prey thrive 
for a long time and then also go extinct. In the second 
route, that can be called (almost) parallel, the prey go 
extinct first, causing a rapid extinction of the predators. 
Assuming large sub-population sizes in the coexistence 
state, we apply a WKB (after Wentzel, Kramers and Bril- 
louin) approximation to the pertinent master equation. 
In this way we evaluate the probabilities of each of the 
two extinction routes and predict the most likely path of 
the sub-populations on the way to extinction. We also 
suggest an effective three-state master equation for the 
evolution of the probabilities to observe the coexistence 
state, the predator-free state and the empty state. 



II. MODEL AND MEAN-FIELD DYNAMICS 

We assume that predators (F - foxes) and prey (R - 
rabbits) are well mixed in space so that spatial degrees of 
freedom are irrelevant. The rabbits reproduce naturally, 
whereas the foxes only reproduce by predation. The foxes 
and rabbits die or leave with constant per-capita rates, 
in general different for each population. We also account, 
in a model form, for the competition of rabbits over re- 
sources, by adding the elementary process 2R — > R that 
becomes more effective at large population sizes. By tak- 
ing into account the competition, this model generalizes 
the classical Lotka-Volterra model (25[. The elementary 
processes and their rates are described in Table HI We 
chose the units of time so that the per-capita death rate 
of the foxes is equal to 1. The large parameter N de- 
termines the typical scaling of the sub-population sizes, 
see below, and we will assume a strong inequality N 3> 1 
throughout the paper. 



Process 


Type of transition 


Rate 


Reproduction of rabbits 


R^2R 


a 


Predation/reproduction of foxes 


F + R^2F 


i/(rw) 


Death of rabbits 


R-^0 


b 


Death of foxes 


F —yO 


1 


Competition among rabbits 


2R -> R 


1/JV 



TABLE I: Predator-prey model 



The mean-field equations for this model are 
R 



(a - b)R - -^rrRF - J-R 2 . 
TN 2N 



F = —RF-F. 
FN 



(1) 



The same model can be reinterpreted to describe the 
spread of an infectious disease in an isolated community. 
Here we re-interpret the rabbits as the susceptibles (S) 
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and the foxes as the infected (I). As in the conventional SI 
model with population turnover [lH [liH HH-tHf . a suscep- 
tible individual can become infected upon contact with 
another infected, while infected individuals are removed 
(recover with immunity, leave or die) with a constant 
per-capita rate. In the conventional SI model the sus- 
ceptibles arrive from outside. In the modified SI model, 
presented here, there are no arrivals from outside, but 
the susceptibles can reproduce by giving birth. They 
are also removed (die or leave), as in the conventional 
SI model. Finally, they compete for resources, 2S — > S, 
so their population size remains bounded. See the ele- 
mentary processes and their rates in Table HU where we 
measure time in the units of per-capita removal rate of 
the infected. The mean-field equations for the modified 
SI model are 



S = (a-b)S - — SI - — S 2 , 
v ' TN 2N ' 

/ = —SI -I. 
FN 



(2) 



These of course coincide with Eqs. (fTJ) upon the change 
of S to R and I to F. 




Process 


Type of transition 


Rate 


Reproduction of susceptibles 


S -» 


2S 


a 


Infection 


I + S 


-> 2S 


1/(TN) 


Removal of susceptibles 


s~ 


¥ 


b 


Removal of infected 


I - 


► 


1 


Competition among susceptibles 


2S - 


-> s 


1/N 



TABLE II: Susceptible-Infected model for an isolated com- 
munity 

For concreteness, we will use the predator-prey nota- 
tion in the following. Introducing the rescaled population 
sizes x — R/N and y — F/N, we can rewrite the mean- 
field equations ([1) as 



x = (a — b)x 



r 



y 



y 



xy 

r 



y- 



(3) 



The mean-field equations are fully characterized by two 
parameters, a — b and T. We will be only interested in 
the case of a — b > and r < = 2(a — 6), when 
Eqs. (|3|) have three fixed points corresponding to non- 
negative population sizes. The fixed point Mi (x\ = 
0, yi = 0) describes an empty system. It is a saddle 
point: attracting in the y-direction (when there are no 
rabbits) , and repelling in the x-direction. The fixed point 
M 2 [&2 = r», ?/2 = 0] describes an established population 
of rabbits in the absence of foxes. It is also a saddle 
point: attracting in the x-direction (when there are no 
foxes), and repelling in a direction corresponding to the 
introduction of a small number of foxes into the system. 
The third fixed point M 3 [x 3 = T, y 3 = T(T* - T)/2] 
is attracting and describes the coexistence state. It is 



FIG. 1: (Color online) The mean field phase trajectories of 
the system for a = 2, 6 = 1 and two different values of T. (a) 
T = 1.8: M3 is a stable node, (b) F = 0.4: M3 is a stable 
focus. 



a stable node for r > Tq 
stable focus for V < Tq. 
a > b. Note also that y 3 
of r. It vanishes at T = 

^2 



= 4(V1 + a - b - 1), and a 
Note that To < for any 
is a non-monotonic function 
and r = r„, and reaches 
a maximum, r^/8, at T = r*/2 corresponding to the 
optimal predation rate. Figure Q] shows two examples 
of mean-field trajectories: for To < T < T» (a) and for 
r < To (b). The characteristic relaxation time scale t r 
of the coexistence state is determined by the real part of 
the eigenvalues of the linear stability matrix at M3. 

Let us compare the mean-field dynamics of this model 
with those of the classical Lotka-Volterra model I2al and 
of the SI model with population turnover [n], [3 113 ■ 
The competition among the rabbits eliminates one un- 
desirable feature of the Lotka-Volterra model: the un- 
limited proliferation of rabbits in the absence of foxes. 
Furthermore, there are no neutral cycles in this model, 
in contrast to the Lotka-Volterra model. The attracting 
fixed point M3, describing a unique coexistence state of 
the rabbits and foxes, appears instead, with either oscil- 
latory (underdamped) or non-oscillatory (overdamped) 
character of phase trajectories approaching it. With ac- 
count of the discreteness of rabbits and foxes, and of the 
stochastic character of their interaction, this leads to an 
exponentially long life-time of the coexistence state, in 
contrast to a much shorter, power-law life-time, charac- 
teristic of the classical Lotka-Volterra model @, [HI ■ 
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On the other hand, the mean-field dynamics ([2j is 
quite similar to that of the conventional SI or SIR model 
with population turnover [TJ, 0, |2f|, except that the 
fixed point Mi (0,0), that is absent in the conventional 
SI model, now appears. As this fixed point is repelling in 
the iE-dircction, its presence docs not make much differ- 
ence in the mean-field description of the coexistence (or, 
in the context of the SI model, of the endemic state of 
the disease in the population). 



III. STOCHASTIC ANALYSIS 



Two extinction routes 



The situation, however, changes considerably when one 
accounts for the stochasticity. Here the empty system 
is an absorbing state describing extinction of both sub- 
populations. Having reached this state via a rare large 
fluctuation, the system will stay there forever. It is this 
empty state that represents the only true steady state of 
this system. (This is different from the conventional SI 
model, where the true steady state describes an estab- 
lished infection- free population.) Importantly, the ab- 
sorbing state at zero can typically be reached, at N ^> 1, 
via one of two routes. The first route can be called se- 
quential: The foxes go extinct first, whereas the rabbits 
thrive for a long time, and then also go extinct. The 
presence of this route is caused by the fact that all fox- 
free states represent absorbing states for the foxes. The 
second route is (almost) parallel: the rabbits go extinct 
first causing a rapid (almost deterministic) extinction of 
the foxes. 

The two different routes to extinction are clearly ob- 
served from Figure [2^ and b that shows two different 
realizations of the stochastic processes, listed in Table 1, 
for the same set of rate coefficients and for the same ini- 
tial conditions. The parameters are the same as in Fig. 
QJ>. In figure a the foxes go extinct first, whereas the 
rabbit population gets established. Then the rabbits will 
also go extinct after a very long time (not shown in the 
figure). In figure b the rabbits go extinct first causing a 
rapid extinction of the foxes. Figure^ and b shows the 
same trajectories in the R, F plane. One can notice in 
Figures [5Jd and [3Jd that the population of rabbits is very 
small close to the time when the foxes go extinct. This 
feature, reproducible in many stochastic simulations, will 
obtain a natural explanation in the WKB theory, see sub- 
section C3. More generally, the WKB theory will enable 
one to compare the probability of each of the two extinc- 
tion routes, to predict the most likely path of the rabbits 
and foxes to extinction, and to evaluate the mean time 
to extinction of each sub-population. 




4000 



FIG. 2: (Color online) Two different stochastic realizations 
of the model for a = 2, b = 1, T = 0.4 and N = 120. 
Shown are the population sizes of rabbits (solid lines) and 
foxes (dashed lines) versus time (measured in Monte Carlo 
steps), (a) The foxes go extinct first, whereas the rabbits get 
established around the fixed point M2 and, after a very long 
time, also go extinct (not shown), (b) The rabbits go extinct 
first, causing a rapid extinction of the foxes. 



B. Master equation and long-time dynamics 

Let P m ,n(t) be the probability to observe, at time t, 
m rabbits and n foxes, where m,n = 0,1,2.... The 
evolution of P m ,n(t) is described by the master equation 

Pm,n — HPm^n — tz[(?71 l)^m-l,)i ^iPm,n] 

+ (l/IW)[(m + l)(n - l)P m+ i,„_i - mnP m ,„] 

+ b[(m + l)P m+ i t „ - raP m ,„] 

+ (n + l)P m , n+ i - nP„ ltn 

+ (l/2N)[(m + l)mP m +i,n - m(m - l)P m ,„],(4) 

where P^j = when any of the indices is negative. There 
are three quantities of interest here describing extinction 
of the sub-populations involved. The probability of ex- 
tinction, by time t, of foxes at any number of rabbits is 
equal to 2ro=o Pmfi{t)- The evolution of this quantity 
in time is described by the equation 



(5) 



m— 



m=0 



directly following from Eq. (Q}. The right-hand side of 
Eq. ([3]) is simply the probability of death of the last re- 
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FIG. 3: (Color online) Same as in Fig. but on the R, F 
plane. 



maining fox at any number of rabbits. 

In its turn, ^2^ =0 Po,n(t) is the probability of extinc- 
tion, by time t, of rabbits at any number of foxes. This 
quantity evolves in time according to 



d 
dt 



:v ^ x x 

E P 0,«W = ^ X>^l,n(*)+&X>,nW- (6) 



n=0 



n=0 



The first and second terms on the right are the prob- 
abilities of predation and of natural death of the last 
remaining rabbit, respectively. Finally, the probability 
Po,o(t) of extinction of both rabbits and foxes by time t 
is described by the equation 



dPo.o 
dt 



bP 1>0 (t) + P oa {t) . 



(7) 



At N 1, and at times much longer than the relaxation 
time t r of the mean-held theory, the quantities that ap- 
pear in Eqs. ©-(0) become sharply peaked around the 
corresponding fixed points of the mean-held theory. The 
structure of these peaks is affected by the presence or 
absence of absorbing states. At long times, the proba- 
bility distribution P m ,o(t) with m > (which describes 
the dynamics of rabbits conditional on prior extinction of 
the foxes) is a one-dimensional distribution peaked at the 
fixed point M 2 of the mean-field theory. The probability 
distribution P m .n(t) with m,n > (which describes the 
long-time dynamics of the coexisting rabbits and foxes) is 
a two-dimensional distribution peaked at the fixed point 



M 3 of the mean-field theory. Finally, the extinction prob- 
ability of both sub-populations Po,o(t) corresponds to a 
Kronecker-delta probability density. Not only the struc- 
ture, but the long-time dynamics of these three distribu- 
tions arc different. To clearly see this point, let us define 
the total "probability contents" of the vicinities of each 
of the fixed points Mi, M 2 and M 3 : 



Vi{t) 



fb,o(i), 

OO 

m— 1 

00 00 



(8) 
(9) 

(10) 



At N > 1 and t > t r the sums in Eqs. © and (JTDJ) are 
mostly contributed to by close vicinities of M 2 and M3, 
respectively. The long-time evolution of Vi, V 2 and V 3 
is described by an effective three-state master equation: 



Vi(t) = r 2 iT 2 (t) + r 31 T 3 (t) , 
V 2 (t) = -r 21 V 2 (t) + r 32 V 3 (t) . 
A(t) = ~(r 3 i+r 32 )V 3 (t), 



(11) 



where r,j is the (yet unknown) transition rate from the 
vicinity of the fixed point i to the vicinity of the fixed 
point j. Let the initial condition correspond to the coex- 
istence state around M 3 : 



[V 1 (0),V 2 (0),V 3 (0)} = (0,0,1). 
Then the solution of Eqs. (fTTj) is 



(12) 



Pi(t) 
Pa(t) 



1 + r 32 e- r ^ t + {r 31 -r 21 )e-^+ r ^ t 



r 2 i - r 3 i - r 32 
r 32 [ e -^+ r ^ t - e- r2lt ] 
r 2 \ - r 3 i - r 32 
_ g— (r3i+r3a)i 



(14) 
(15) 



Once the transition rates r 3 i, r 32 and r 2 \ are known, 
Eqs. (fT3|) -(fT5 ]) provide a valuable "coarse-grained" de- 
scription of the stochastic system in terms of the long- 
time evolution of the probabilities to observe the coexis- 
tence state around M 3 , the fox- free state around M 2 and 
the extinction state at M\. The coexistence probabil- 
ity V 3 (€) goes down to zero exponentially in time. The 
probability Vi(t) of extinction of both sub-populations 
increases monotonically with time, exhibiting two differ- 
ent exponents, and ultimately reaches 1. Finally, the 
probability of the fox- free state V 2 (t) first increases with 
time, reaches a maximum and then goes down to zero. At 
intermediate times, i r « t <C min [l/r 2 i, 1/(^31 + ^32)], 
Eqs. ([15 ]) - ([IB ]) read 



Vi{t) a r 31 t, V 2 (t)~r 32 t, 
V 3 (t) ~ I - {r 31 + r 32 )t. 



(16) 



During this stage of the dynamics, Vi(t) and V 2 (t) grow 
linearly with time, whereas the transition 2 — > 1 does 
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not show up. What happens at longer times? The tran- 
sition rates are usually widely different in magnitude. 
According to our WKB calculations, presented below, 
r 2 i < r 31 + r 32 . Then, at times t > l/(r 3 i + r 32 ), 
Eqs. ([I3 ]) -([T5 |) become 



v 2 (t) 



1 



r 32 e 



»~3i + r 32 



0. 



<"32 



(17) 

(18) 
(19) 



2_R, R — >• and 2i? — > i?, where the predators are absent 
(l8| . The eigenvalue r 3 i + r 32 is also exponentially small 
in N, see below. It corresponds to the quasi-stationary 
distribution 7r mi „ at n > 0, sharply peaked at the fixed 
point M 3 . Our main effort in the following will be to 
evaluate, with exponential accuracy, the transition rates 
r 3 i and r32 that contribute to the eigenvalue r 3 i + r32. 
As all the effective transition rates r\j are exponentially 
small in N S> 1, we can drop the right-hand side in the 
eigenvalue problem 



Hltr, 



(r 3 i + 7-32)^77 



n > 



(23) 



Now there is a probability current from state 2 to state 1, 
describing extinction of rabbits in the absence of foxes, 
but the transition rates r 3 i and r 3 2 are still present in 
the equations. 

It is common, see e.g. Ref. [2j, to characterize the 
extinction risk of a stochastic population in terms of its 
mean time to extinction. From Eqs. (|13p ~ (|15p . the mean 
time to extinction of foxes is 



and arrive at the quasi-stationary equation 



Hn ri 



, m > 0, n > . 



C. WKB approximation 

1. General 



(24) 



T F = I dtt[P!(t) + P 2 (t)} = dttP 3 (t) 
'0 Jo 

1 



r 3 i + r 32 



(20) 



whereas the mean time to extinction of both sub- 
populations is 

t rf = / dttA(t) - 7 32 r . (21) 
Jo r 2 i(r 3 i +r 32 ) 

How is this dynamics encoded in the spectral proper- 
ties of the the linear operator H, see Eq. (g])? Let Aj be 

the eigenvalues and ir m \n the eigenstates of H. These are 
defined by the relations 



Hir (t} 



-AV l) 



so that 



-Pm,n (^) 



1,2,.. 



where the constants C% are determined by the initial con- 
dition P TOn (0). Let us order the eigenvalues so that 

Ai < A 2 < A 3 < The true (empty) steady state 

corresponds to the only zero eigenvalue: Ai = 0, whereas 
7Tm,n = S m odnO- Equations dT3l) - (fT5)) and the inequality 
?"2i < ^31 +^32 imply that, at N 3> 1, the next two eigen- 
values are A2 = r 2 i and A3 = r 3 i + r 32 . The quantity r 2 i 
was found in Refs. [13, [HI, and it is exponentially small 
in N: 




exp 



-27V 



b + b In ■ 



(22) 

The corresponding eigenvector 7r m! o^nO is the quasi- 
stationary distribution of the one-population model R — > 



For N 3> 1, Eq. ([M]) can be approximately solved via 
the WKB ansatz 



TTm.n = eXp[-JVS(x, J/)] , 



(25) 



where S is assumed to be a smooth function of the contin- 
uous variables x = m/N and y = n/N . The use of WKB 
approximation for finding stationary or quasi-stationary 
solutions of master equations with a discrete state space 
was pioneered in Refs. [29T - [32l | . By now this approxi- 
mation, in the space of population sizes or in the mo- 
mentum space, has become a standard tool in the anal- 
ysis of rare lar ge f luctuations of stochast ic p opulations 

SHE [Ull M M M MM IS- Plugging 
Eq. (J2S) into Eq. (JH) and Taylor expanding S to first 
order around (2, y), we arrive at a zero-energy Hamilton- 
Jacobi equation H(x,y,d x S,d y S) = 0, with the Hamil- 
tonian 

H(x, y,p x ,Py) = ax(e Px — 1) + bx(e~ Px — 1) 

2 

+ ^( e Py-P* _ 1) + y ( e -P« _ 1) + £_ ( e -P. _ 1). (26) 

The trajectories are given by the Hamilton's equations 
for the "coordinates" x and y and conjugate "momenta" 
p x and p y : 



bx + f_ ) e -p. _ ^ e P«-P, 

2 / r 



Pa 



(b + x)(l - e- p *) + I (1 - e p «"^) + o(l - e p *), 

(27) 



+_(l-e P »^), 



where we are only interested in zero-energy trajectories. 
The (zero-energy) invariant hyperplane p x = p y = cor- 
responds to the mean- field dynamics ([3]). The invariant 
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hyperplanes x = and y = correspond to the rabbit- 
free and fox-free dynamics, respectively. The Hamilto- 
nian problem is characterized by three independent pa- 
rameters a, b and T. 

The Hamiltonian flow ([27]) has five zero-energy fixed 
points: 

Mi = (0,0,0,0), 

m 2 = (r,, 0,0,0), 

m 3 = [r,r(r*-r)/2,o,o], 

Fx = [0,0,ln(6/a),0], 

f 2 = [r*,o,o,-in(ryr)]. (28) 



The one-dimensional quasi-stationary distribution 
7i"m,o^n,o, corresponding to the long-lived population of 
rabbits in the absence of foxes, corresponds to another 
instanton-like trajectory: the one connecting the fixed 
points Mi and Fi 0, [3 ■ The accumulated action 



S21 



Pxdx 



M-2 



b + bhx- 



(30) 



yields the transition rate r 2 i which agrees, u p to a pre- 
exponential factor, with Eq. (J22J. See Refs. [U, [3 for 
detail. 



The zero-momentum fixed points Mi , M 2 and M3 corre- 
spond to the three fixed points of the mean field equa- 
tions ([3]), so we keep the same notation for them. The 
two other fixed points, F\ and F 2 , are fluctuational fixed 
points describing a fox-free state at a non-zero number 
of rabbits (F 2 ) and an empty system {Fx). Fluctuational 
fixed points have a non-zero p x or p y component and 
appear in a broad class of stochastic population mod- 
els exhibiting extinction in the absence of an Allee ef- 
fect 0, M ED, m M HE US Hi, M M Hi- They 

play an important role in the calculations of the quasi- 
stationary distributions and the mean time to extinction. 
It is mostly their presence that makes the WKB theory of 
population extinction distinct from the WKB theory of 
noise-induced switches between different states that are 
stable in the mean-field theory [39j. 

To determine S(x,y) in Eq. (|2"5|) , one should calcu- 
late the action accumulated along the (zero-energy) ac- 
tivation trajectory in the phase space of the Hamiltons 
equations ([27]) that exits the fixed point M3 and ends at 
(x,y)- 

S(x,y) = / p x dx+p y dy, 
Jm 3 

and then minimize the result with respect to p x and p y at 
the end point (x,y). To evaluate the transition rates r 3 i 
and r3 2 , we need to evaluate ir m ,n at points [x = 0, y = 0) 
and [x — T*,y = 0), respectively. As in many other 
stochastic population models, there are no phase trajec- 
tories that would start at M3 and end at fixed points 
Mi or M 2 . There are, however, instanton-like activation 
trajectories that start, at t — —00, at M3 and enter, at 
t = 00, the fixed point Fi or F 2 , respectively. The accu- 
mulated actions S31 and S32 along these instantons, 



S31 



p x dx+pydy and ^32= / p x dx+p y dy, 

M 3 



yield, with exponential accuracy, the transitions rates r 3 \ 
and r 32 : 



r-31 

T32 



cxp(-NS 31 ) 
exp(-NS 32 ) 



2. Analytic result for S32 near bifurcation T = T» 

Because of the lack of an independent integral of mo- 
tion in addition to the Hamiltonian, the canonical equa- 
tions (|2"T|) are non-integrable analytically for a general set 
of parameters a, b and T. Close to the bifurcation T = r„, 
however, the fixed points M3 and F 2 become very close 
to each other, and the action 5*3 2 can be calculated per- 
turbatively by exploiting the time separation property 
uncovered in Refs. [Hi EH- Let us define the bifurcation 
parameter 5=1 — r/r*. For <5 <§; 1, the fixed points 
M 3 and F 2 become M 3 ~ [^(1-^,^/2,0,0] and 
F 2 ~ (r*,0, 0,— S), Motivated by these expressions, we 
assume the following scalings: x = T*+6qi, y = STlq 2 /2, 
p x = S 2 pi, and p y = Sp2, where qi, q 2 , p\ and p 2 are 
(non-canonical) rescaled variables. The equations of mo- 
tion become, in the leading order in 5, 



y (gi + r*g 2 ), 



1i 

q_2 = <%2 1 1 - 
Pi 

P2 



r* 



r* 

Y (pi - 92P2) , 



= 6 [pi 



q\P2 



1 



P2 - P 2 



(31) 



The dynamics of q\ and p\ is fast compared to that of 
q 2 and p 2 . As a result, qi and p\ adiabatically follow 
the slowly varying q 2 and p 2 : qi = — T*q 2 , p\ = q 2 p 2 . 
Plugging these expressions in the equations for g 2 and 
P2, we obtain equations 

qi = Sq 2 (1 + 2p 2 - q 2 ) , P2 = $P2 (2<7 2 - 1 - P2) , (32) 

derivable from the reduced "universal" Hamiltonian 

#r(<?2,P 2 ) = Sq 2P 2(p2 -92 + 1) [IS M M H ■ This 

problem is easily soluble, and one obtains 



S 2 T 2 



(29) 



(q 2 -l)dq 2 



(r*-r) : 



(33) 
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FIG. 4: The numerically found instantons M3F1 and M3F2 
for a = 2, b — 1 and T = 1.8, when the fixed point M3 
is a node. Shown are the icy-projections (a) and the p x p y 
projections (b). 



3. Numerical calculations of S31 and S32 

For a general set of parameters a, b and T, one can 
find the instantons M 3 Fi and M 3 F 2 numerically: either 
by shooting [Tol . [Tl| , or by iterating the equations for x 
and y forward in time, and equations for p x and p y back- 
ward in time [U, H3, 5(J ■ Here we employed the shooting 
method and explored different parts of the parameter 
space a, b and T. Figure 2] shows typical examples of 
the numerically found instantons M3F1 and M3F2 in the 
case when M3 is a node. Figures [5] and [6] refer to the 
case when M3 is a focus. As one can see, the instanton 
M3F2 is qualitatively similar, in both regimes, to that 
found in the conventional SI model [ll|. The instanton 
M3F1 is different, as it corresponds to extinction of both 
sub-populations, absent in the conventional SI model. 

Note that Fig. [6^l resembles Fig. [3k,, and Fig. [5^, resem- 
bles Fig. In particular, the instanton M3F2 passes 
close to the point (0,0). This explains the feature ob- 
served, for the same set of parameters a, b and T, in 
stochastic simulations, see Fig. As in other stochastic 
systems exhibiting rare large fluctuations, one should ex- 
pect that by averaging over a a large number of stochastic 
trajectories, conditional on a given extinction route, one 
will obtain the corresponding instanton up to an error 
that vanishes asJV^oo. 

Figured shows the accumulated actions S31 and S32 
for fixed a — 2 and b = 1 and different values of T (param- 
eterized by (5 = 1 — r/r«). The first observation is that 

531 > S'32 for all 8, implying that the effective transition 
rate r^2 is exponentially greater than r^\. As expected, 

532 vanishes at 8 — > and 8 —¥ 1 and has a maximum 
close to 8 = 1/2, or V = T+/2. This behavior follows that 
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FIG. 5: The numerically found instanton M3F1 for a = 
2, b — 1 and T — 0.4, when the fixed point M3 is a fo- 
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cus. Shown are the xj/-projection (a) and the p x p y projection 
(b). 



8 




0.2 0.4 0.6 0.8 1 

5 



* 


(b) 







0.2 0.4 0.6 0.8 1 

d 



FIG. 7: (Color online) Accumulated actions S31 (solid line) 
and S32 (dashed line) as functions of 8 = 1 — T/T* for a fixed 
6=1 and two different values of a: 2 (a) and 3 (b). The 
dotted lines show the analytical asymptote from Eq. (|33[1 . 
The diamonds denote the values of S21 from Eq. (f30)l . 



of the mean-field population size of foxes in the coexis- 
tence state. S31 behaves differently: it has a maximum 
at 8 (that is, at r — > r„), goes down monotonically 
with an increase of <5, and tends to zero as S approaches 
1 (that is, r approaches 0). This behavior can be qual- 
itatively understood from the T-dependence of the fixed 
points, see Eq. ([25)) . There is, however, one surprising 
feature here. As 8 approaches zero, the fixed point M3 
becomes closer and closer to the fixed point M2. One 
might expect, therefore, that the instanton M3F1 and 
the accumulated action S31 would approach the instan- 
ton MiF\ and the action 5*2 1, respectively. This is not 
what happens. Extrapolating the numerically calculated 
values of S31 toward 8 = 0, we obtain 831(8 = 0) ~ 0.48. 
This is considerably less than 621 (8 = 0) = 0.6137..., 
obtained from Eq. (|30l) . see Fig. [7^,. Furthermore, the nu- 
merically found instantons M3F1 at small 8 are markedly 
different from the instanton A/2-F1 for which y = p y = 0. 
Not only the instantons M3F1 exhibit relatively large in- 
termediate values of y and p y , but the maximum value 
of p y along the instanton actually increases as 8 goes 
down. That is, as 8 — > 0, or T — > r», the fluctuations in 
the number of foxes play an important role in the joint 
extinction of the rabbits and foxes. 

We also observed these features for other values of pa- 
rameters a, b and T from the coexistence region < 
r < T*. For example, Fig. [7b shows the ^-dependence 
of the accumulated actions S31 and S32 for a = 3 and 
6=1. Here the two actions are very close to each 
other in a broader region of <5, but the double inequality 
S32 < S31 < S21 still holds. Furthermore, it continues 



to hold when a is only slightly greater than b (not shown 
here). 

When amplified by the very large factor N 3> 1 in 
the exponent, see Eq. (|29p . the double inequality S32 < 
S31 < S21 leads to 

rai < r 3 i < r 32 , (34) 

implying that the sequential extinction route (foxes first, 
rabbits second) is much more likely than the (almost) 
parallel extinction route. Using the inequality r 3 2 3> rai, 
we can further simplify Eqs. (fT7]) - ([2Tj) . In particular, the 
mean time to extinction of foxes becomes simply Tp ~ 
I/V32, whereas the mean time to extinction of both sub- 
population is trf ~ l/r2i. Here the quantities Tp and 
trf are determined by the "kinetic bottlenecks" of the 
effective transitions, as expected. 

However, for a sufficiently high predation rate, and not 
too large N, the transition rates r^x and r 3 2 are compara- 
ble, and Eqs. (fl7]) - (PT|) should be used in their complete 
form. 



4- Proximity of instantons M3F1 and M3F2 at small F 

As is evident from Fig. the actions S31 and S32 ap- 
proach each other closely for sufficiently large 8, that is 
for high predation rates. The reason for it becomes clear 
upon comparison of the instantons M3F1 and M3F2 for 
sufficiently small T, see Fig. [SJ One can see that the in- 
stanton M3F2 almost coincides with the instanton M3F1 
until, close to the fixed point Fi, it abruptly changes its 
direction and goes toward the fixed point F2 . On the lat- 
ter segment of the trajectory the values of y and p x stay 
close to zero, so the contribution of this segment to the 
action S32 is negligible. In other words, the most likely 
route to extinction of the foxes in this regime of parame- 
ters involves a drastic decrease in the number of rabbits: 
almost all the way to their extinction! Then, a few re- 
maining rabbits start reproducing, at a very small num- 
ber of foxes, and the rabbit population gets reestablished. 
Note that, in the WKB formalism, the reestablishment 
of the rabbits, accompanied by extinction of the foxes, 
occurs via the "fluctuational" fixed point F2 , rather than 
via directly approaching the mean-field point M2. 



5. Two modifications of the model 

One central result of this work is the strong inequal- 
ity P31 *C ^32 observed in most of the parameter space 
that we explored. That is, the sequential extinction route 
(first the predators then, after a long time, the prey) is 
usually much more likely than the (almost) parallel route. 
In other words, the predators are usually much more vul- 
nerable extinction-wise than their prey. In retrospect, 
this feature is hardly surprising: As the predators are 
dependent on their prey for survival, they can be at best 
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FIG. 8: (Color online) xy (a) and p x p y (b) projections of the 
instantons M3F1 (solid lines) and M3F2 (dashed lines) for 
a = 3, b = 1 and T = 1.5. For these parameters S = 0.625, 
and the actions S31 and 5*32 are very close, see Fig. [7b. 



as resilient extinction-wise as the prey. (This also gives 
an intuitive explanation to the proximity of the instan- 
tons M3F2 and M3F1 at high predation rates, reported in 
the previous subsection.) Will this feature hold for other 
predator-prey models? To get insight, we introduced two 
minor modifications of our model, each of them endowing 
the predators with more resilience but still keeping them 
dependent on the prey for survival. In one modification 
we added a direct reproduction process F — > 2F with 
rate d < 1. This can model a large amount of additional 
small prey for the foxes, say mice. The direct reproduc- 
tion reduces the effective death rate of the foxes. Still, 
as long as d < 1, the foxes go extinct deterministically 
in the absence of rabbits. In another modification of the 
model we replaced the predation process F + R —> 2F 
by a more efficient one F + R — > 3^. Our numerics 
showed that, for both of these modified models, the in- 
equality 53.1 > 53,2 continues to hold. Furthermore, in 
both models we observed, for sufficiently high predation 
rates, the convergence of the instantons M3F2 and M3F1 
to each other until a close vicinity of the fixed point F\ is 
reached. A difference between the original model and the 



modified ones is that, as we improve the conditions for 
the predators, this convergence occurs at lower predation 
rates. 



IV. CONCLUSIONS 

We considered a simple stochastic predator-prey model 
in its coexistence region and observed that this model ex- 
hibits two different routes of extinction of each of the pop- 
ulations: the sequential and the (almost) parallel. This 
multiplicity of the extinction routes can be conveniently 
accounted for by an effective three-state master equation 
for the probabilities to observe the coexistence state, the 
predator-free state and the empty state. The WKB ap- 
proximation yields the effective transition rates between 
these three states, as well as the most likely paths of 
the sub-populations to extinction: the instantons. We 
showed numerically that the parallel extinction route is 
usually much less likely than the sequential one, implying 
a great robustness of the prey against predation. For a 
sufficiently high predation rate, however, the two routes 
may have comparable probabilities. 

The rest of parameters being fixed, our model predicts 
an optimal predation rate so that the mean time to ex- 
tinction of the predators is maximum. Not surprisingly, 
this optimal predation rate is close to that for which the 
quasi-stationary population size of the predators is maxi- 
mum. A surprising result is that, for sufficiently high pre- 
dation rates, the optimal path to extinction of the preda- 
tors almost coincides with the optimal path to the joint 
extinction of the predators and prey until a point (corre- 
sponding to an almost zero size of each sub-population) 
is reached where the two optimal paths depart from each 
other: the predator population continues moving toward 
extinction whereas the prey population reestablishes it- 
self. That is, for a high predation rate, the predators are 
more likely to reach extinction by consuming nearly all 
the prey. This phenomenon appears to be robust and in- 
dependent of details of the predator-prey model, as long 
as predators still need prey for their survival. 

Finally, all our results can be reformulated in terms 
of the SI epidemic model for an isolated community, see 
Table 2, where both the infected, and the susceptible 
populations are prone to extinction. 
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